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Multilayer infrastructure is often interdependent, with nodes in one layer depending on nearby 
nodes in another layer to function. The links in each layer are often of limited length, due to the 
construction cost of longer links. Here, we model such systems as a multiplex network composed of 
two or more layers, each with links of characteristic geographic length, embedded in 2-dimensional 
space. This is equivalent to a system of interdependent spatially embedded networks in two dimen¬ 
sions in which the connectivity links are constrained in length but varied while the length of the 
dependency links is always zero. We find two distinct percolation transition behaviors depending on 
the characteristic length, £, of the links. When £ is longer than a certain critical value, £ c , abrupt, 
first-order transitions take place, while for C < Cc the transition is continuous. We show that, though 
in single-layer networks increasing £ decreases the percolation threshold p c , in multiplex networks 
it has the opposite effect: increasing p c to a maximum at £ = £ c . By providing a more realistic 
topological model for spatially embedded interdependent and multiplex networks and highlighting 
its similarities to lattice-based models, we provide a new direction for more detailed future studies. 


I. INTRODUCTION 

Several models have been proposed for spatially em¬ 
bedded networks nHiDi- In lattice-based models, links 
are only formed to nearest or next nearest neighbors. In 
random geometric models, links are formed to all neigh¬ 
bors within some distance HUE]. In models of power 
grid topology, links are formed with the m nearest neigh¬ 
bors, statically m or as a generative model M- Some 
models utilize a cost function [T51 - H8] or a characteristic 
distance distribution □ms to determine link lengths. 
The model which we study here has spatiality expressed 
via characteristic link lengths. We utilize exponentially 
distributed link lengths, similar to the Waxman model 

HD- 

Interdependent networks have been studied mainly on 
random topologies where analytic calculations are pos¬ 
sible [24l - l29l . However, since many real-world com¬ 
plex systems are embedded in space, it is important 
to understand the properties of interdependent networks 
with topology reflecting the dimensionality of the space 
[I6tl30| . This is particularly important when dealing with 
critical infrastructure which is heavily influenced by spa¬ 
tial constraints El HU- An important first step in 
that direction is the model presented by Li et al. [33] 
which models the networks as lattices and includes de¬ 
pendency links which are of limited geographic length 
(described by the parameter r). This model was shown 
to include a number of surprising properties including 
three regimes of phase transitions depending on the value 
of r. For r < r c , the phase transition is second-order 
and appears to be in the same universality class as 2- 
dimensional lattices pAj. The percolation threshold in¬ 
creases as r increases and reaches a maximum value at r c , 
where the transition switches to first-order EH [3Sj. For 
r c < r < oo the transition is first-order and characterized 
by a spreading process, withp c decreasing monotonically. 
For r = oo, the transition is an abrupt simultaneous first 
and second order transition mm- The model was stud- 



FIG. 1: Examples of real-world networks with links of 
characteristic length. We examine the distribution of the 
geographic lengths of the edges in both the European power 
grid [22] (1851 edges) and the inter-station local railway lines 
in Japan EH (20745 edges). These networks have links of 
characteristic length and longer links are exponentially un¬ 
likely, as indicated by the linear drop on the semi-logarithmic 
plot. To compare the two datasets, we rescale the lengths so 
that they are measured in units of their own minimum length, 
which we determine as the peak of the distribution (mode 
length). The normalization value (l = 1) corresponds to 3.7 
km (power) and 1.0 km (rail). The characteristic length as 
determined by the mean is 4.8 km (power) and 1.2 km (rail). 
As the slope of the exponential fit, it is 3.3 km (power) and 
2.0 km (rail). The Japan local railway data is formed from 
the complete railway network from [23] with bullet train lines 
and internal station tracks removed. 

ied under partial dependency with r = oo and it was 
shown to have first-order transitions for any fraction of 
dependency EH]- When there is partial dependency and 
finite dependency link lengths, r c is shown to increase 
as the fraction of dependent nodes decreases, diverging 
at q = 0 i3B., j35]. This model was also studied for gen- 
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FIG. 2: Spatially embedded multiplex networks, (a) The nodes occupy regular locations in two-dimensional space while 
the links in each layer (blue and green) have lengths that are exponentially distributed with characteristic length £ = 3 and 
are connected at random, (b) Cascading failures in multiplex networks. In the first stage the mutual giant connected 
component (MGCC) consists of the entire network since all of the nodes are in the giant component of both layers (blue links 
and green links). An initial attack on Node 3 causes it and its links to fail. This detaches Node 5 from the giant component of 
the green links, and in the next step it and its links fail. After the failure of Node 5, Node 6 is no longer in the giant component 
of the blue links. After Node 6’s failure we find that the remaining nodes are in the giant component of both layers. We note 
that the MGCC is not simply the intersection of the giant components in the separate layers. For example, using such an 
approach one would conclude that Node 6 is in the MGCC after the failure of Node 3, which is not the case. 


eral networks formed of interdependent lattices [?D] for 
interdependent resistor networks with process-based de¬ 
pendency [35] and in the presence of healing [4Tj. 

However, in many real-world systems, the length of the 
dependency links may not be longer than the length of 
the connectivity links. For instance, in the example of 
the power grid and communications network, it is un¬ 
likely that a communications station will skip the nearer 
power stations and depend on a power station that is far¬ 
ther away. Also, real-world networks including the power 
grid do not have uniform link lengths like a lattice, but 
rather have links of a characteristic length, consistent 
with an exponential distribution (see example in Fig. [IJ 
[30] . To address these two issues, here we model interde¬ 
pendent networks where every node has a bi-directional 
dependency link with the nearest geographic node in the 
other network. We treat each pair of nodes as single 
nodes in a multiplex network, where the links in each 
layer are different but of the same characteristic length 
We thus interpret each node as a geographic entity 
(e.g. a city or neighborhood) which is linked via two 
types of links (e.g. electricity and communications) to 
other nodes. Each node requires both of its constituents 
to function and each constituent requires connectivity 
within its layer. 

Our main focus in this paper is to examine the role 
of the characteristic length (£) of the links on the ro¬ 


bustness of the multiplex. We find that, though increas¬ 
ing £ decreases p c in single networks-making them more 
robust-it has the opposite effect on multiplex networks. 
Increasing £ increases p c for multiplex networks until a 
critical length £ c where p c is maximal. At £( c , the per¬ 
colation transition changes to first-order and the multi¬ 
plex network is susceptible to spreading cascading fail¬ 
ures similar to the ones observed in interdependent lat¬ 
tices [3U [311 [35J 32] • 

By demonstrating that comparable critical behavior 
emerges in more realistic topologies, we show that the 
critical behavior demonstrated in previous lattice-based 
models is not limited to the specific implementation or 
the lattice topology but is rather a generic property of 
interdependent spatially embedded networks. Further¬ 
more, by providing a topological model that more closely 
matches real-world systems, we provide a more realistic 
topological framework for future studies of interdepen¬ 
dent critical infrastructure. 


II. MODEL 

To model connectivity links of characteristic length in 
each layer, we construct the network as follows. We be¬ 
gin by assigning each node an ( x , y) coordinate with in¬ 
tegers x, y £ [0,1,..., L). To construct the links in each 
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(a) single-layer (b) multiplex 

FIG. 3: Percolation of (a) single-layer and (b) multiplex networks with links of characteristic geographic 
length. The fraction of the network in the largest connected component (Poo) as a function of the fraction of nodes remaining 
after random removal of a fraction 1 — p nodes, (a) In single networks, the transition is always second-order and the critical 
behavior above p c is the same for all finite values of £ (( < L). The value of p c decreases quickly and monotonically with £ 
(£ > 1). (b) In multiplex networks, the transition is comparable to single networks for £ = 0.2 ( p c ta 0.5927) but p c increases 
as £ increases. This continues until the maximal value is reached (at C , c ) and the transition becomes first-order. The case of 
£ = 1000 is very close to purely random and has p c ~ 2.4554/(fc) = 0.61385. Each line represents an individual realization for 
(k) = 4 and L = 4000. 


layer, we select a source node at random with coordinates 
(xo,Vo) and draw a length l according to P(/) ~ e~ 1 ^. 
We choose the permitted link length (dx,dy) which is 
closest to fulfilling l = y/dx 2 + dy 2 , select one of the 
eight length-preserving permutations (dx -f-)- — dx 7 dy -f-> 
— dy,dx -O- dy) uniformly at random and make a link to 
node (xi,j/i) with X\ = xq + dx , y\ = yo + dy. This 
process is executed independently in each layer and is 
continued until the desired number of links (N(k)/ 2) is 
obtained. For simplicity, we use the same characteristic 
length £ and average degree (k) in each layer. However, 
because they are constructed independently, the links in 
each layer are different (as demonstrated in Fig. 2a and 
Fig. § and this disorder enables the critical behavior 
which we describe below. 


We then perform site percolation by removing a frac¬ 
tion 1 — p of the nodes from the system and finding the 
mutual giant component |24| . When a node is removed, it 
causes damage to nodes in both layers due to the connec¬ 
tivity links in each layer, which are also removed. How¬ 
ever, since the connectivity links are not the same in both 
layers, there will be nodes that are connected to the giant 
component in one layer but are disconnected in the other 
layer. Since the node functionality requires connectivity 
in both layers, such nodes will fail, causing further dam¬ 
age in the system. This leads to the cascading failures as 
demonstrated in Figs. |2b| and Fig. [5 which are similar 
to those described in (24, 31113 Ell SMI- 


III. RESULTS 


Since we are not aware of a discussion of the perco¬ 
lation properties of this topology for single networks, 
we briefly describe those properties here and in the Ap¬ 
pendix. In the limit of £ —> 0, the only permitted links 
will be to nearest neighbors (because links of length < 1 
are not accessible) and a square lattice is recovered. As 
such, in the case of (k) = 4, we recover the standard 
2-dimensional percolation behavior with p c ss 0.5927 
(Fig. 3a). When slightly longer (i.e., next nearest 


neighbor) links are allowed, the system becomes slightly 
less robust as discussed in the Appendix. As ( increases 
further, the robustness increases and in the limit £ —► oo, 
all lengths are equally likely to be drawn and the topol¬ 
ogy reverts to purely random (Erdos-Renyi topology) 
with p c = 1 /(k) = 0.25 (see Fig. |3a[) Thus, similar to 
rewiring probability in the original small world model [1] , 
we have a single parameter £ which allows us to smoothly 
transition from lattice to random topology. For all values 
of <C, a single network undergoes a second-order transition 
(Fig. 4c and 4d). 


In multiplex networks, where connectivity to the giant 
component in both layers is required, cascading failures 
emerge Ei nans HU- For C ~ 0, large cascading failures 
do not emerge. This is because the multi-layer structure 
is mostly redundant and the difference between connec¬ 
tivity in one layer or both layers is negligible. However, 
once ( becomes long enough (above ( c ), intensive cascad¬ 
ing failures emerge and the system undergoes an abrupt, 
first-order transition, similar to the transition in inter¬ 
dependent lattices (33] > as shown in Fig. [4] and [5j As ( 
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FIG. 4: The effect of the characteristic link length ( on the percolation threshold and size of the giant 

component at criticality, (a-b) The percolation threshold in multiplex networks increases until it reaches a peak at and 
then decreases slightly to its asymptotic value of 2.4554/(fc). (c-d) The size of the order parameter at p c . Single-layer 
networks always have an order parameter close to zero at p c , which indicates a second-order transition. Multiplex networks 
have a second-order transition for low values of £ but at Q this jumps to a large fraction of the system size, indicating a 
first-order transition. All panels have L = 4000 and are plotted based on at least 10 realizations of each system. The averages 
and raw data are plotted for all systems. 


becomes even longer, p c decreases and slowly approaches 
its asymptotic value of 2.4554/(fc) as known from inter¬ 
dependent Erdos-Renyi networks [23], (Fig. 4a and 4b). 

In interdependent lattices the mechanism of the first- 
order transition for dependency links of large finite length 
is a propagating spinodal interface [31 05} [321 02 ■ On 
a microscopic level, the first-order transition takes place 
by the emergence of a hole due to random fluctuations 
of characteristic length £ (the percolation correlation 
length) which then propagates through the system. This 
propagation is enabled by the cascade dynamics (cf. Fig. 


2b) and dependency links which relay the damage caused 


by the hole into a concentrated area around the hole’s 
edge. It would seem that in order for this phenomenon to 
take place, long dependency links are needed and longer 
connectivity link lengths would be insufficient to sustain 
damage propagation. The reason that the dependency 
links have a stronger influence on robustness is that in 
order for a node to fail due to connectivity, all of its links 
to the giant component need to fail whereas with depen¬ 
dency, an otherwise well-connected node will fail if the 
single node that it depends on fails. 

We find that, in fact, this is not the case. Connectiv¬ 
ity link lengths which are above a critical length, £ c , but 
much smaller than the total system size are sufficient to 
cause the cascading failures observed in lattice models, 
even when the dependency links have zero length. In¬ 
deed, the first order transition lacks scaling behavior just 
above p c (Fig. |3b[) a nd is characterized by a slow spread¬ 
ing process (Fig. [5]), just like interdependent lattices de¬ 
scribed in Refs. [31I23I3H1- However, the maximum p c 
in this system is lower than the comparable lattice sys¬ 
tem: « 0.64 ((fc) = 4 here) vs. « 0.74 (interdependent 
lattices). This is because, in order for the connectivity 
links to effectively relay the damage from the hole, the 
system must be closer to criticality. Indeed, only when 


the average degree of each node is 1 do the connectivity 
and dependency links have the same effect upon it. 

Surprisingly, the highly localized topology which 
emerges from spatial embedding makes the system more 
susceptible to cascading failures. This is due to the fact 
that, because the damage from an emergent (or induced 
[421 ) hole is relayed by the cascade dynamics to the neigh¬ 
borhood of its interface, the nodes near the edge of the 
hole are far more likely to become disconnected. This 
is related to work on information diffusion in social net¬ 
works, where it was shown that high modularity makes 
viral cascades more likely to occur due to the increased 
likelihood of multiple exposure to the information @3- 

[51] . 

In single layer networks, p c decreases monotonically 
as ( increases from ( « 0.5 to the limit of £ = oo. In 
contrast, in multiplex networks, p c increases until (j c and 
then decreases monotonically thereafter. The peak in 
p c (0 is due to the fact that the size of the critical hole 
that is needed to trigger the transition scales linearly with 
( [32] and the size of emergent holes above the percola¬ 
tion threshold, £(p), decreases with p [33]. This would 
indicate that the smaller £ is, the smaller the critical 
hole needs to be and that p c would increase monotoni¬ 
cally as £ decreases. However, when ( < ( c there is not 
enough space between the emergent hole and the extent 
of the damage propagation (£) for the network to dis¬ 
integrate and the small emergent holes remain in place 

[31ES1M1S2- 

For interdependent lattices with dependency links of 
finite length, the single network case is recovered when 
the dependency links have length zero (r = 0) [331 • 
spatially embedded multiplex networks the same limit 
is recovered as ( —>■ 0 due to overlapping links. The 
fraction of common connectivity links between two in¬ 
terdependent networks or between two layers in a multi- 
























5 



(a) t = 30 




(b) t = 100 



(c) t = 200 (d) t = 300 


FIG. 5: Dynamic evolution of cascading failures. Here 
we show the nodes of the multiplex network, colored black 
if functional and white if not. After the emergence of a 
large enough hole (due to random fluctuations), the dam¬ 
age is relayed predominantly to the area around the edge of 
the hole (due to the finite link lengths) which leads to the 
nodes around the interface becoming disconnected and prop¬ 
agation of the damage until the system disintegrates. Similar 
dynamics have been observed in interdependent lattices with 
dependency lengths of finite length [Ml EH- In this figure 
L = 4000, C = 15, (k) =4. 



FIG. 6: The fraction of overlapping nodes. The fraction 
of overlapping nodes is determined as the number of common 
links across both layers divided by the total number of links 
in each layer. When ( « 0, the overlap is maximal and the 
networks are identical (for (k) = 4, as in this figure). As £ 
increases, the fraction decreases with an exponent of ss —1.8, 
see text for discussion. 


the deviation from the continuum calculation is due to 
the fact that with the discretization of space that we 
introduce, links that would otherwise have been distinct 
are unified and the critical exponent is reduced from —2 
to ~ —1.8. 

Unlike studies of random multilayer networks with 
overlap 03J 00] [S2; ( 155] . decreased overlap alone is not 
enough to enable the critical behavior observed here. It 
is only the combination of the disorder (as indicated by 
decreased overlap) with the spatially embedded links that 
enables the distinctive first-order transition which we ob¬ 
serve here. 


plex network is called intersimilarity [43j 04] or overlap 
[521 (53]. The cascading failures and abrupt transitions 
which characterize interdependent networks decrease as 
overlap increases. In the limit of total overlap, they dis¬ 
appear altogether because the system is composed of two 
exact copies of the same network and the cascade shown 
in Fig. [2b] does not take place at all. 

In multiplex networks with links of characteristic 
length, the extent of the overlap can be estimating by 
considering the probability that, given the same source 
node, two links lead to the same target node. In the con¬ 
tinuum limit this is proportionate to the probability that 
the links have the same length and the same direction. 
The system is isotropic by construction so the directional 
condition is simply l/27rr, the size of a ring of radius r. 
We obtain 


P(overlap) 


P 2 {r) 

2nr 


dr 


1 f°° e~ 2r ^ 
C 2 A 2?rr 


1 


(1) 

We find that the scaling in our system is scale-free, with 
an exponent of ss —1.8 (Fig. [6|. We hypothesize that 


IV. DISCUSSION 


Previous models of spatially embedded interdependent 
networks [50] S3 I581BU1 021 have used two-dimensional 
lattices with dependency links connecting nodes from 
one network to the other. The dependency links were 
also affected by the spatial embedding via the restriction 
that they have length of up to r, a system parameter. 
This model led to many important results, but left sev¬ 
eral important issues unaddressed. First, the topology of 
real-world spatially embedded networks is not lattice-like 
or even strictly planar and it was not clear that results 
derived on lattices would accurately describe real-world 
topologies. Second, the assumption that dependency 
links are longer than connectivity links does not corre¬ 
spond with what we would expect from critical infras¬ 
tructure: It is not reasonable to expect a communications 
station to get power from a distant power plant and not 
the one nearest to it. Here we address these problems by 
modeling spatially embedded interdependent networks as 
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(a) (k) = 3 (b) (fc) = 4 

FIG. 7: Dependence of p c on £ for single networks The 

percolation threshold p c drops quickly as a function of £ and 
by £ ss 10 it is already very close to p c = 1 /(k), the value 
from Erdos-Renyi networks. 

multiplex networks where the dependency relationship is 
to the nearest node in the other layer and the connectiv¬ 
ity links are of finite characteristic length but not uniform 
or regular. 

We find that the most important features of the lattice- 
based models are reproduced by our new model: second- 
order and first-order transitions depending on the link 
length, a first-order transition characterized by the emer¬ 
gence and spreading of a hole and substantially higher 
vulnerability compared to single-layer networks. Fur¬ 
thermore, p c has a maximal value at the £ value (£ c ) that 
separates the two types of transitions. This validates pre¬ 
vious work based on lattices and also shows a new way 
forward for the modelling of critical infrastructure and 
other spatially embedded multilayer networks. 
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Appendix: Percolation threshold in single networks 
and in multiplex and single networks for £ < 3 

In single networks, for £ > 1 , p c decreases monotoni- 
cally with increasing (. This is evident in Fig. [TJ where 
p c is shown to rapidly approach 1 /(fc), the Erdos-Renyi 
value, as £ is increased Thus, we find that the 

effect of increasing £ on single networks is to make them 
more robust, the opposite of its effect in multiplex net¬ 
works. The behavior of single and multiplex spatially 
embedded networks is different for 0 < ( < 3. Since, this 
is not the main focus of this study and since it does not 
affect the critical phenomena, we have avoided discussing 


(a) ( k > = 3 (b) (k) = 4 

FIG. 8: Percolation behavior for very short character¬ 
istic lengths. The percolation threshold increases slightly 
for all systems with very small £. For single networks it 
rapidly falls again but for multiplex networks, p c remains 
higher for longer before falling to its £ = 0 level. This is be¬ 
cause in single networks, the only effect is the relative weak¬ 
ness of lattices with complex neighborhoods. In multiplex 
networks, this effect is compounded with the weakness en¬ 
gendered by clustering [56} [K7] and overlap j43[ |44} [52} 020 • 


it in the main text. When £ = 0, only the shortest pos¬ 
sible links will be drawn (those with l = 1). Thus, the 
topology of the network layer will be a pure lattice for 
(k) = 4 and a diluted lattice for ( k} < 4. However, for 
0 < C < 3, some of the nearest neighbor links are ex¬ 
changed for next nearest neighbor (diagonal) links which 
are less robust than nearest neighbor only links. This 
is a well-known result from bond percolation on lattices 
with complex neighborhoods, where it was found that 
(k c ) for a next-nearest neighbor lattice (z = 8) is higher 
than (k c ) for a nearest neighbor lattice (z = 4) [53]. For 
instance, for £ ss 0.6 the network layer is a mixture of 
approximately 85% nearest neighbor links and 15% next 
nearest neighbor links. As £ increases, the links are less 
redundant and tend to strengthen the network, as seen 
in Fig. [7] and Fig. [8] for single networks. 

The robustness is more severely impacted in the mul¬ 
tiplex case. This is due to the effects of clustering and 
disorder. Combining nearest and next-nearest neighbor 
links leads to high clustering (number of triangles) which, 
though having only a minor impact on the robustness 
of single networks, substantially weakens interdependent 
networks HJ E7|. Furthermore, the increased disorder 
decreases the fraction of overlapping links substantially 
(~ 50% as in Fig. [6]). These effects combine in the 
multiplex case to cause a substantial increase in p c for 
0.3 < £ < 1- At ( « 1, the links are no longer redun¬ 
dant, clustering decreases, and the system becomes more 
stable again. However, it is important to note that even 
though p c is increased in this interval, the transition is 
still second-order. The first-order transition only takes 
place when the links are substantially longer (£ > £ c ), 
and the spreading process described above can take place. 
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